#include "dyn_quad_qRomberg.h"

qRomberg::qRomberg(Functor &_func,
                   double _a,
                   double _b,
                   double _tol,
                   int _ptsInterpolation,
                   int _maxIterations)
    : Quad(_func,_a,_b,_tol)
{
    setPtsInterpolation(_ptsInterpolation);
    setMaxIterations(_maxIterations);
}

double qRomberg::doQuadrature()
{
    // We prepare vectors to hold the successive
    // refinements and their respective stepsizes.
    colvec s(maxIterations), h(maxIterations+1);
    h(0) = 1.0;                                     // The first stepsize is the unit.

    PolyInterpolation polint(h,s,ptsInterpolation); // We prepare the polynomial interpolation.
    TrapzDriver t(func,a,b);                          // And the quadrature driver.

    for (int i=1;i<=maxIterations;i++)
    {
        // Compute the next refinement.
        s(i-1) = t.computeNextRefinement();

        // After a few steps, we can use interpolation to
        // enhance convergence.
        if (i>ptsInterpolation)
        {
            // We compute the interpolation value.
            polint.setXValues(h), polint.setYValues(s);
            double sInterpolated = polint.rawInterpolation(i-ptsInterpolation,0.0);

            // We check if we have converged.
            if (fabs(polint.getError()) <= tol*fabs(sInterpolated))
                return sInterpolated;
        }
        h(i) = 0.25*h(i-1);
    }
    throw("Failed quadrature.");
}
